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ABSTRACT 

A flexible maximum-entropy component separation algorithm is presented that accom- 
modates anisotropic noise, incomplete sky-coverage and uncertainties in the spectral 
parameters of foregrounds. The capabilities of the method are determined by first 
applying it to simulated spherical microwave data sets emulating the COBE-DMR, 
COBE-DIRBE and Haslam surveys. Using these simulations we find that is very diffi- 
cult to determine unambiguously the spectral parameters of the galactic components 
for this data set due to their high level of noise. Nevertheless, we show that is possible 
to find a robust CMB reconstruction, especially at the high galactic latitude. The 
method is then applied to these real data sets to obtain reconstructions of the CMB 
component and galactic foreground emission over the whole sky. The best reconstruc- 
tions are found for values of the spectral parameters: Td — 19 K, ad — 2, — —0.19 
and /3syn = —0.8. The CMB map has been recovered with an estimated statistical 
error of ~ 22/zK on an angular scale of 7 degrees outside the galactic cut whereas the 
low galactic latitude region presents contamination from the foreground emissions. 

Key words: methods: data analysis-cosmic microwave background. 



1 INTRODUCTION 

Observation of the cosmic microwave background (CMB) 
radiation is one of the most powerful tools of modern cos- 
mology. A handful of experiments - Boomerang (Netterfield 
et al. 2002), MAXIMA (Hanany et al. 2000), DASI (Halver- 
son et al. 2002), VSA (Grainge et al. 2003), CBI (Mason 
et al. 2003), ACBAR (Kuo et al. 2004), Archeops (Benoit 
et al. 2003), WMAP (Bennett et al. 2003a) - have aheady 
reported the measurement of CMB fluctuations at subde- 
gree scales, allowing tight constraints to be placed on cosmo- 
logical parameters. Moreover, current and future CMB ex- 
periments will measure the fluctuations with unprecedented 
resolution, sensitivity, sky and frequency coverage. Most no- 
tably, these include the WMAP mission by NASA (that will 
continue to take data in the next years) and the Planck mis- 
sion by ESA (to be launched in 2007), both of which will 
provide all-sky multifrequency observations of the CMB. 

When measuring the microwave sky, however, one does 
not only receive the cosmological signal but also emis- 
sion from our own Galaxy, thermal and kinetic Sunyaev- 
Zeldovich emissions from clusters of galaxies and emission 
from extragalactic point sources. In addition, instrumental 



noise and possibly some systematic effect will contaminate 
the data. Therefore, our capacity to extract all the valuable 
information encoded in the CMB anisotropics will depend 
critically on our ability to separate the cosmological signal 
from the other microwave components. 

Different methods have been proposed in the literature 
to perform this component separation. Some techniques at- 
tempt to separate and reconstruct all the components at 
the same time, such as Wiener filtering (Bouchet, Gispert 
& Puget 1996, Tegmark & Efstathiou 1996), the maximum- 
entropy method (Hobson et al. 1998, Stolyarov et al 2002) 
or blind source separation (Baccigalupi et al. 2000, Maino 
et al. 2002, Delabrouille et al. 2003). Another approach is 
to extract only the component of interest from the data, as 
is the case, for example, when one is searching for emission 
from extragalactic point sources or the Sunyaev-Zel'dovich 
effect in clusters. Methods for performing such a separation 
include the use of the mexican hat wavelet (Cayon et al. 
2000, Vielva et al. 2001a, Vielva et al. 2003) matched fllters 
(e.g. Tegmark & de Oliveira-Costa 1998), scale-adaptive fil- 
ters (Sanz, Herranz & Martinez-Gonzalez 2001, Herranz et 
al. 2002a,b,c) - see also Barreiro et al. (2003) for a compar- 
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ison of the performance of these filters -, the McClean algo- 
rithm (Hobson & McLachlan 2003), the Bayesian approach 
proposed by Diego et al. (2002) or the blind EM algorithm 
of Martinez-Gonzalez et al. (2003). 

In the present work, wo will focus on the maximum- 
entropy method (MEM) for reconstruction of all compo- 
nents simultaneously. This technique has been successfully 
applied in reconstructing microwave components from sim- 
ulated Planck data in a small patch of the sky (Hobson et 
al. 1998, hereafter H98). This work was extended to deal 
with point sources (Hobson et al. 1999) and also combined 
with the mexican hat wavelet (MHW; Vielva et al. 2001b). 
Recently, the MEM and MEM+MHW algorithms have been 
adapted to deal with spherical data at full Planck resolution 
by Stolyarov et al. (2002) (hereafter S02) and Stolyarov et 
al. (2003) respectively. The algorithm is capable of analysing 
the vast amount of data expected from the Planck mission 
by working in harmonic space, assuming independence be- 
tween difi'erent harmonic coefficients. Although this assump- 
tion is not necessarily true, the method is very successful in 
performing a full-sky component separation. Nevertheless, 
it would be desirable to find a method that can perform 
this task without the assumption of independent harmonic 
modes, and can be applied directly in the space where the 
data have been observed. This would also allow one to intro- 
duce the properties of the noise in a more straightforward 
manner and to deal simply with incomplete and arbitrary 
shape sky coverage. 

Another shortcoming of the standard MEM approach to 
component separation is that the spectral dependence of the 
microwave components needs to be known a priori. Although 
this is the case for the CMB and the kinetic and thermal SZ 
effects, the spectral behaviour of the galactic foregrounds is 
only approximately known. Moreover, it varies with position 
and frequency. The point sources, in particular, can cause 
problems, since each source will have a different spectral 
behaviour. This last point has been solved by combining 
MEM with the MHW, which is optimised for the detection 
of point sources. Nevertheless, the standard approach still 
lacks a way to estimate the spectral dependence of the diffuse 
components present in the data. 

In this paper, we present a maximum-entropy compo- 
nent separation method that works in both real and har- 
monic space and is able to deal with many of the prob- 
lems mentioned above. Our analysis also includes a thorough 
study of the properties of the reconstructions to estimate the 
(average) spectral parameters for the galactic components. 
Unfortunately, the price one has to pay for this flexible MEM 
is to make the method much slower than the harmonic MEM 
used in S02. To illustrate the performance of tfic algorithm, 
we apply it first to simulated spherical data and then to a set 
of real spherical data available during the development of the 
algorithm, including the COBE-DMR, COBE-DIRBE and 
Haslam maps, from which we reconstruct the CMB emission 
and galactic foreground components. 

The outline of the paper is as follows. In Section 2 
we outline the flexible maximum-entropy component sep- 
aration method. The spherical microwave data set that we 
have analysed is described in Section 3. Sections 4 and 5 
are devoted to the analyses of simulated (including an ar- 
tificially low noise case) and real data respectively. Finally, 
we present our discussion and conclusions in Section 6. 



2 THE SEPARATION ALGORITHM 

In this section we outline the basics of our separation 
method, focusing particularly on the differences between the 
algorithm used in this work and traditional harmonic-based 
MEM component separation. For a more detailed deriva- 
tion of MEM, see Hobson et al. (1998) and Stolyarov et al. 
(2002). 

2.1 The problem 

Our aim is to reconstruct the CMB amisotropies and the fore- 
ground components in the presence of instrumental noise 
from multifrequency microwave observations. We will as- 
sume that point sources have been previously subtracted 
using the MHW, or other filtering technique, and that their 
residual contribution is negligible. If we observe the mi- 
crowave sky in a given direction x at n/ frequencies, we ob- 
tain a n/-dimensional data vector that contains the observed 
temperature fiuctuations in this direction at each observing 
frequency, plus instrumental noise. The observed data at the 
z/th frequency in the direction x can be written as 

Np no 

j=l p=l 

where Np denotes the number of pixels in each map and Uc 
is the number of physical components to be separated. As 
is usual for the MEM algorithm, we make the assumption 
that each of the components can be factorised into a spatial 
template (sp) at a reference frequency vo and a frequency 
dependence encoded in F^p. The function B^, accounts for 
the instrumental beam and e,^ (x) corresponds to the instru- 
mental noise at frequency v and position x. 

2.2 Harmonic-space MEM 

If the beam is circularly symmetric and du is measured over 
the whole sky, it is convenient to write the former equation 
in harmonic space as 

ric 

dem = ^ RV^'im + ''tm, (2) 
p=l 

where we have adopted the usual notation for spherical har- 
monic coefficients fim = J^^ dQYi^{yi) f {x) in which Yim{x) 
is a standard spherical harmonic function. Therefore, 
and correspond to the spherical harmonic coefficients 
of the data and the noise at the i^th frequency respec- 
tively, whereas are the harmonic coefficients of the p 
component at the reference frequency. The response matrix 
R'^^ = BiFvp (where Bi are the harmonic coefficients of the 
z/th observing beam) determines the contribution of each 
physical component to the data. Using matrix notation, we 
have for each mode 

dlrn = R^a^m + ^im, (3) 

where d^™,, a.em and etm are column vectors of dimension 
71/, ric and rif complex components respectively, whereas 
the response matrix contains Uf x ric elements. 

If one neglects any correlations between different {£, m) 
modes, the reconstruction can be performed mode-by-mode. 
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which greatly simphfies the problem. This corresponds to 
assuming that the emission from each physical component 
and the instrumental noise are isotropic random fields on the 
sky. Although, in reality, this is not the case, and so correla- 
tions do exist between difltcront modes, the modo-by-modc 
harmonic-based MEM produces excellent reconstructions. In 
this approximation, the a priori covariamce structure of the 
different components is assumed to take the form 

{aimS-e'm') = CiSteiSmrn' , (4) 

where we have made the additional simplifying assumption 
that the covariance matrix does depend only on £. Analo- 
gously, the (cross) power spectra for the noise can be written 
as 

{eimee'm') = ^eSuiSmm' (5) 

Thus, although the algorithm assumes statistical isotropy, it 
can straightforwardly include any a priori knowledge of the 
(cross) power spectra of the physical components and the 
noise at different observing frequencies. 

As explained in H98, any a priori power spectrum in- 
formation can be provided to the MEM algorithm by intro- 
ducing a 'hidden' vector him that is a priori uncorrelated 
(see H98) and relates to the signal through 

a-im = L^h^m, (6) 

where the lower triangular matrix is obtained by perform- 
ing the Cholesky decomposition of the cross power spectra 
Ci of the physical components at uq. Therefore, the sepa- 
ration problem can be solved in terms of the hidden vector 
hern and the corresponding a.em are subsequently found us- 
ing (6). Note that L itself can be itcrativcly determined by 
the MEM (H98). That is, one can use an initial guess for L 
to obtain a first reconstruction and subsequently compute 
the power spectra of those reconstructions, which axe used 
as a starting point for the next iteration, until convergence 
is achieved. 

As shown in S02, harmonic-based MEM finds the best 
reconstruction for the sky by minimising mode-by-mode the 
function (for a detailed derivation see S02): 

${hem) = x^(h<m) - aS{hem, m), (7) 

where is the standard misfit statistic in harmonic space 
given by 

X^(h^m) = {dim - 'Re'Lehem)^'N'['^ {dem. - R^L^hf^) (8) 

and S'(h^m,in) is the cross entropy (the form of which is 

given in H98) of the complex vector h£,„ and the model 
m, to which him defaults in absence of data. The regular- 
ising parameter a can be estimated in a Bayesian manner 
by treating it as another parameter in the hypothesis space 
(see §2.3.2). 

Working in harmonic space and neglecting correlations 
between different modes vastly reduces the computational 
requirements of the component separation problem . Instead 
of performing a single minimisation of ~ 2ncimax parame- 
ters, one performs ~ £max minimisations of 2nc parameters, 
which is much faster. In addition, working in harmonic space 
provides us with a simple manner of introducing (cross) 
power spectra information in the algorithm through the L 
matrix. 



2.3 Flexible MEM 

Although the advantages of working in harmonic space are 
clear, there are also some shortcomings. In addition to 
the necessity of neglecting coupling between different {£, m) 
modes, a full and regular coverage of the sky is needed in or- 
der to keep the harmonic transformations simple. The prop- 
erties of the noise are also somewhat diluted in harmonic 
space and it is not obvious how anisotropic or correlated 
noise is affecting the algorithm. Finally, it is not straightfor- 
ward how to combine data taken in different spaces (such as 
ID scans, intcrfcromctric data, incomplete spherical data, 
etc.). Therefore, it would be desirable to develop a MEM 
algorithm that is able to deal with all these issues, but still 
keeps as many advantages of the harmonic-based MEM as 
possible. Unfortunately, this is a non-trivial issue. 

As a first approach, we have implemented a MEM al- 
gorithm that combines the space where the data have been 
taken and harmonic space. The best sky reconstruction is 
obtained by performing a single minimisation, with respect 
to the whole set of parameters h = {hf^}, of the function 

VP(h) = xI(h)-aS(h,m), (9) 

where Xd is evaluated in the space where the data have been 
taken. Note that, in this case, there is no need to trans- 
form the data into harmonic space. Moreover, the properties 
of the noise are well defined in data space and anisotropic 
noise can be easily included in the analysis. This is simply 
achieved by including the relevant noise dispersions in the 
noise covariance matrix of the misfit statistics Xd- Also, in- 
complete sky coverage or galactic masks can be taken into 
account, since the Xd can be calculated summing only over 
a portion of the data. This method also allows the combi- 
nation of different set of data by combining their Xd val- 
ues. The cross entropy, however, is still calculated in har- 
monic space by summing over the entropy S(hf,„,m) for 
each mode, thereby preserving the straightforward intro- 
duction of power spectrum information into the algorithm. 
The reconstruction is thus still obtained in harmonic space, 
making necessary an inverse transform from the harmonic 
modes to real space (which is usually much simpler than 
the forward transform) . The form of Xd > £^s a function of the 
hidden variables htm, for spherical pixelised data is given in 
Appendix A. Since we are not performing a mode-by-mode 
minimisation, correlations between modes may be taken into 
account by the algorithm. 

Unfortunately, the price that one has to pay for such a 
flexible method is the need to perform a single minimisation 
of order ^ ricimax parameters, which makes the method 
many times slower than harmonic MEM. Even so, we believe 
that it is worth exploring the possibilities of this algorithm, 
since it could be extremely useful in many applications that 
can not currently be performed by harmonic MEM. 



2.3.1 Newton-Raphson minimisation 

In order to perform the minimisation of the function given 
in (9), we use a Newton-Raphson (NR) iterative algorithm. 
In the full NR method, the sky modes at iteration i + 1 are 
obtained from their values at iteration i via 

h'+i=h^-7(ff)-V, (10) 
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where 'loop gain' 7 is a parameter of order unity, whose 
optimal value is determined through the method itself, g' 
is an n — vector (where n is the number of parameters to 
be minimised) containing the gradient of the function with 
respect to h evaluated at h = h\ and H is the curvature 
(or Hessian) matrix of dimension nxn evaluated at h* . The 
form of the first and second derivatives of S are given in 
H98 and those of for spherical pixelised data are given in 
Appendix A. 

Unfortunately, with this scheme, one needs to calculate 
and invert the Hessian matrix with dimension nxn, which 
is not feasible even for low-resolution data. Nevertheless, the 
main part played by the Hessian matrix in the NR method 
is to provide a scale length for the minimisation algorithm. 
Thus, the success of the algorithm does not require a very 
accurate determination of the Hessian (indeed, many min- 
imisation algorithms achieve very good results using only 
gradient information). In practice, therefore, we calculate 
only an approximation to the Hessian matrix in the NR al- 
gorithm. We have found that good results are obtained by 
approximating the Hessian by a block diagonal matrix with 
blocks of size ric xuc, and setting the remaining elements to 
zero. This corresponds to assuming that 

unless {l,m) = {£',m') = {£",m") and k = k' = k" , where 
the indices k refers to the real or imaginary part of him and 
c to each of the considered components to be reconstructed. 
In this case, the inversion of the matrix can be performed 
block-by-block, which simplifies the problem enormously. 
Note that the value of each hem. mode still depends on the 
values of all the rest of the modes through the gradient vec- 
tor (see Appendix A). 



2.3.2 Determining a 

Another important issue is the determination of the regular- 
isation parameter a. As explained in H98, a can be obtained 
in a fully Bayesian manner by including it as another param- 
eter in the hypothesis space, and maximising the Bayesian 
evidence with respect to it. In particular, one finds that it 
must satisfy 

aS'(h) =n-QTr(M-^), (11) 

where M = G^^/^HG"^/^ and G is the (diagonal) metric 
of the image space (see H98 for details). Note that in order 
to solve this implicit equation we also need to operate with 
the Hessian matrix. To mate this task feasible, we again 
approximate H hy a block diagonal matrix with blocks ric x 
He- This allows one to find a (nearly) optimal value for the 
Q parameter. 



2.3.3 Error estimation 

Following HQS and S02, we can also obtain an estimation of 
the covariancc matrix of the reconstruction errors Saim on 
the harmonic modes a^^, as well as the dispersion of the 
residual map for each component. In particular, by again 
approximating by a block diagonal matrix, we find 



{Sa.tmSa.im) = ((»£,„ — atm) {aim — aim) ) = LiHf.^Ll (12) 

where Hem is the ric x ric block of the Hessian matrix corre- 
sponding to the (£, m) mode evaluated at hem- As shown by 
S02, the former equation allows us to obtain the residuals 
power spectrum at component p simply by 

e 

m=-e 

where (|5a^^|^) is given by the pth. diagonal entry of the 
diagonal matrix L^H^l|. Finally, we can estimate the dis- 
persion eest of each residual map as: 

e 

2.4 Estimating spectral behaviour 

In order to apply the MEM algorithm, one needs to assume 
that the spectral behaviour of all the components to be re- 
constructed is known and spatially constant. This is the case 
for the CMB and the SZ (flVx-ts, but it is not true for the 
diffuse galactic components. Although there has been con- 
siderable effort in recent years to study the galactic fore- 
grounds at microwave frequencies, there are still uncertain- 
ties in the knowledge of the frequency dependence of free- 
free, synchrotron and dust emissions (as well as in their spar 
tial distribution). In addition, this frequency dependence is 
expected to vary across the sky. To study the synchrotron 
emission, the most extensively used data has boon the 408 
MHz map of Haslam et al. (1982). More recently, Reich & 
Reich (1986) and Jonas et al. (1998) measured the emission 
of the northern (at 1420 MHz) and southern (at 2326 MHz) 
galactic hemispheres respectively. At these frequencies the 
synchrotron emission dominates the total galactic emission, 
and therefore these data sets are very useful in characteris- 
ing this foreground. The combination of these maps allows 
one to estimate the frequency dependence of the synchrotron 
in this frequency range. Reich & Reich (1988) obtained an 
average value of the synchrotron spectral index of -2.7 (as- 
suming a power law in temperature units) in the northern 
hemisphere. Giardino et al. (2002) generated an spectral in- 
dex map for the whole sky, whose average value is also - 
2.7. More recently, the WMAP team (Bennett et al. 2003b) 
found that the synchrotron power law in the WMAP range 
frequency is relatively flat inside the galactic plane with an 
index of -2.5 whereas it steepens to —3 outside the Galaxy. 

Thermal dust emission is usually modelled by a grey- 
body whose emissivity depends on the physical properties 
of the materials that constitute the dust (see e.g. Desert, 
Boulanger & Puget, 1990; Banday & Wolfendale, 1991). 
Recently, there has been an effort to produce a dust tem- 
plate, as well as to determine its frequency dependence at 
microwave frequencies. In particular, Schlegel, Finkbeiner 
& Davis (1998) produced a dust map at 3000 GHz using 
IRAS and COBE-DIRBE data. They also characterised the 
dust emission in the 1250-3000 GHz range with a grey body 
law with an emissivity of Qd = 2 and a varying spatial 
temperature with values around 17-21 K. More recently, 
Finkbeiner, Davis & Schlegel (1999) proposed an improved 
two-component dust model using the IRAS, COBE-DIRBE 
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and COBE-FIRAS data. The so-called cold component is 
characterised by a mean temperature of 9.4 K and an emis- 
sivity of 1.67 whereas the spectral parameters of the hot 
component take values of 16.2 K and 2.7 respectively. 

The least well known of the galactic foregrounds is the 
free-free emission. Usually, it is estimated through two differ- 
ent tracers: the thermal dust emission and the Ha emission. 
Up to a few years ago, no surveys were available and ther- 
mal dust was normally used to produce free-free templates, 
taking into account the correlation found for instance be- 
tween the COBE-DIRBE and COBE-DMR data (Kogut et 
al. 1996). Recently, several Ha surveys have been produced: 
the Virginia Tech Spectral line Survey (VTSS) of the north- 
ern hemisphere (Dennison, Simonetti & Topasna 1998), the 
Wisconsin H- Alpha Mapper (WHAM) that covers a large 
fraction of the sky (Reynolds, Haffner & Madsen 2002) and 
the Southern H-Alpha Sky Survey Atlas (SHASSA) of the 
southern hemisphere (Gaustad et al. 2001). Even more re- 
cently, all-sky free-free templates have been generated com- 
bining the information of the different surveys (Dickinson, 
Davies & Davis 2003, Finkbeiner 2003). The spectral de- 
pendence of the free-free emission is normally modelled as 
a power law (in temperature) with spectral index around 
-2.16 (e.g. Kogut et al. 1996, Smoot 1998). 

There are also uncertainties in the number of com- 
ponents that contribute to the microwave sky. In fact, an 
anomalous galactic emission at low frequency, which is well 
correlated with the thermal dust one, has been found by sev- 
eral authors (de Oliveira-Costa et al. 1997, Leitch et al. 1997 
and Kogut 1999). A possible candidate for this component 
has been proposed by Draine & Lazarian (1998): electric 
dipole emission coming from rapidly rotating dust grains 
("spinning dust"). The first statistical evidence for spinning 
dust was given by de Oliveira-Costa et al. (1999). Finkbeiner 
et al. (2002) found two tentative detections of this emission 
in two (out of ten) small areas of the sky. However, Bennett 
et al. (2003b) established that the spinning dust contribu- 
tion to the WMAP frequencies is clearly subdominant. 

In summary, in order to apply our MEM algorithm, and 
taking into account all the uncertainties in the knowledge of 
the galactic foregrounds, we need a way to model their spec- 
tral behaviour, which makes use of our prior knowledge of 
the foregrounds as well as of the data themselves. Addition- 
ally, it would be desirable to be able to accommodate spatial 
variations of the spectral parameters. 

One possibility to determine the (spatially constant) 
frequency parameters from the data themselves would be 
to use an iterative approach, assuming a known spectral 
law for the components. The procedure would be as follows: 
(i) reconstruct the microwave components with MEM using 
an initial guess for the unknown spectral parameters; (ii) 
use the reconstructions as templates to fit for the spectral 
parameters and the normalisation of the components (by 
minimising x^); (iii) use these new parameters as a starting 
point to run MEM again; (iv) repeat the procedure until 
convergence. We have tested this method on our set of sim- 
ulated data, but unfortunately it did not always converge to 
the correct values. However, these data have a low signal- 
to-noise ratio and we need to fit for several spectral pa- 
rameters at the same time, which leads to degeneracies. In 
addition, MEM will drive the reconstructions towards the 
templates that best fit the initial (incorrect) spectral pa- 



rameters, which may be far from the templates that fit the 
correct parameters. All these factors make it very difficult to 
estimate the frequency dependence of the components with 
this method for our set of data. Nevertheless, if high qual- 
ity data arc available, or if a single component needs to be 
fitted, this method should be further investigated. 

A different approach is to run MEM for different sets of 
spectral parameters and try to infer from the reconstructed 
components which parameters give the best results and are, 
therefore, closest to the truth. One could naively look at 
the estimated errors and pick the reconstructions with the 
lowest values for these errors. However, the estimation of 
the errors depend on the chosen spectral parameters. Basi- 
cally, it gives the statistical error of the reconstruction but it 
does not take into account uncertainties in the values of the 
spectral parameters. Therefore, if our guess of the frequency 
dependence is incorrect, the error estimation of the recon- 
structions is not reliable. Jones et al. (2000) showed that the 
dust frequency dependence could be estimated from Planck 
simulated data of small patches of the sky by minimising 
the of the reconstructions. They assumed a (spatially 
constant) grey body law for the dust emission with two un- 
known parameters (the dust temperature and emissivity), 
which they were able to fit from the data by looking at the 
minimum of the x of the reconstructions. However, varying 
the spectral index /3 of the synchrotron or free- free emissions 
(assuming a power law I oc ly^) had little effect on the value 
of the and could not be determined in this way. This was 
due to the fact that both the synchrotron and free-free emis- 
sions have a low amplitude with respect to the rest of the 
components in the considered region of the sky and therefore 
the data could not provide enough information to fit these 
two components. In this case, the reconstructions of these 
two emissions were lost, but those of the other components 
were little affected. This application shows that, depending 
on the characteristics of the data (sky coverage, resolution, 
signal-to-noise ratio, frequency coverage, etc.) we may not 
be able to estimate some of the spectral parameters using 
just the x^ value. Therefore, in order to estimate all the spec- 
tral parameters we may need to use also information coming 
from other variables, such as the entropy or the cross cor- 
relations between the reconstructed CMB and the galactic 
components. We have used such an approach to determine 
empirically the best reconstructions, which is explained in 
detail in §4.2.2. 



3 SPHERICAL MICROWAVE DATASETS 

The only spherical CMB dataset available during the de- 
velopment of the flexible MEM algorithm was COBE-DMR 
data. In order to enable the reconstruction of different com- 
ponents of emission, we also make use of the COBE-DIRBE 
and Haslam maps. Each of these maps exist in the HEALPix 
pixelization (Corski et al. 1999) with A'^sidc = 32, which cor- 
responds to 12288 pixels of size 107 arcminutes. Since each 
of these datasets has been described in detail elsewhere, we 
just give a brief description in this section. The data maps 
in M.Jy/sr are shown in Fig. 1. 

The COBE-DMR data consists of three frequency 
maps: 31.5, 53 and 90 GHz (each of them obtained by the 
combination of two channels) with a resolution of ~ 7 de- 
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Figure 1. Data maps used in our analysis in units of MJy/sr. The left column shows the three COBE-DMR channels (from top to 
bottom): 31.5, 53 and 90 GHz. The top right panel gives the Haslam map and the middle and bottom right panels correspond to two 
of the COBE-DIRBE channels (1249 and 2141 GHz). The Haslam and COBE-DIRBE maps have been plotted in a non-linear scale to 
allow the visualization of the structure of the maps outside the galactic centre 



grees and a signal-to-noise ratio of around 2 per 10 degree 
patch. The COBE-DMR data have an anisotropic noise pat- 
tern that can be easily taken into account with our method. 
As an illustration we show the noise dispersion per pixel for 
the 53 GHz channel in Fig 2. The COBE beam is well char- 
acterized by the curve given in Fig. 3 (Wright et al. 1994). 
For comparison, a Gaussian beam with 7 degree of full width 
half maximum (FWHM) is also shown. 

The COBE-DIRBE experiment observed the brightness 
of the full sky at ten wavelengths from 1.2 to 240 microns as 
well as mapping linear polarization at 1.2, 2.2 and 3.5 mi- 
crons. In our analyses below we use the data at the two low- 
est frequencies (1249 and 2141 GHz). The maps have been 
degraded down to Nside ~ 32 and smoothed with a Gaussian 
beam of FWHM 2.4 times the pixel size (i.e., FWHM =263.8 
arcmin). The COBE-DIRBE noise is also anisotropic, but 




Figure 2. Dispersion noise per pixel of the COBE-DMR, 53 GHz 
map in MJy/sr. 
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Figure 3. Window function of the COBE-DMR instrument (solid 
line). For comparison a Gaussian beam of FWHM=7° is also 
shown (dotted line). 

since we have significantly degraded the resolution of the 
data, its level is very small. The COBE-DIRBE beam has a 
resolution of 0.7 degrees and can be approximately modelled 
by a top-hat. At these frequencies, (thermal) dust emission 
dominates and therefore these data help the algorithm to 
extract the dust component. 

The Haslam map gives the emission of the whole sky at 
408 MHz and has been obtained by combining four different 
surveys (Haslam ot al. 1982). It has an effective resolution of 
51 arcminutes. Recently, Finkbeiner, Davis & Schlegel (pri- 
vate communication) have reprocessed the original Haslam 
map to provide a destriped, point-source subtracted map, 
which we have used for our analysis. We have degraded its 
resolution down to the HEALPix resolution of Ngide = 32 
and smoothed it with a Gaussian beam of FWHM=263.8 ar- 
cmin. At this frequency the emission is dominated by galac- 
tic synchrotron, providing MEM with excellent information 
to trace this component. 



4 ANALYSIS OF SIMULATED DATA 

Before applying of our flexible MEM component separation 
method to real data, we have checked its performance on 
simulated datasets. In order to simulate our sets of data 
(COBE-DMR, COBE-DIRBE and Haslam maps) we have 
assumed that the sky contains, CMB, frcc-froo, synchrotron 
and thermal dust emissions. The four simulated maps at the 
reference frequency of 50 GHz are given in Fig. 4 in units of 
IjK (thermodynamical temperature). 

4.1 The simulations 

The CMB realization has been produced with the help of the 
HEALPix package and has been constrained to have a power 
spectra compatible to the one derived from the COBE-DMR 
data (Bennett et al. 1996) at low £'s and from Archeops 
(Bcnoit ct al. 2003) for higher multipoles. The synchrotron 
emission has been modelled by a power law 7syn oc 
and the template at the reference frequency has been ob- 
tained by extrapolating the Haslam map with this law. The 
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Table 1. Assumed observational parameters for the simulated 

data 





Haslam 




DMR 




DIRBE 


Preq (GHz) 


0.408 


31.5 


53 


90 


1249 2141 


Resolution (°) 


4.5 


~ 7 


~ 7 


~ 7 


4.4 4.4 


Au/u 


0.04 


0.07 


0.04 


0.05 


0.28 0.40 



dust emission has been obtained by extrapolating the IRAS- 
DIRBE map* of Schlegel, Finkbeiner & Davis (1998) using 
a single grey body component with dust cmissivity of 2 and 
temperature of 18 K. Finally, the free-free component has 
been modelled with a power law of Ig oc z/~°^^. As template 
for the free-free we have used the WHAM survey (Reynolds, 
Haffncr & Madscri 2002), which provides an He map of most 
of the sky. The empty part of the sky has been filled in by 
copying on it another region of the survey. The map has 
been normalized to have a dispersion of 185 fiK at 50 GHz 
and smoothed with a 7 degree Gaussian beam. This is sev- 
eral times higher than expected using estimations from the 
Ha emission (e.g. Finkbeiner 2003), but it is intended to 
account for the excess of emission found in the data whose 
origin is uncertain. 

The separate components have been combined in order 
to simulate data according to the characteristics given in 
Table 1. For the COBE-DMR data we have used the true 
beam shape given in Fig. 3 and added Gaussian pixel noise 
according with the anisotropic pattern of the maps. Since we 
aim to reconstruct only four components (and in particular 
a single dust component) we have simulated only the low- 
est frequency COBE-DIRBE channel. We have ignored the 
beam of this experiment since it is small in relation to the 
pixel (and the subsequent smoothing). We have added the 
corresponding anisotropic noise at the COBE-DIRBE map. 
Signal plus noise have then been smoothed with a Gaussian 
beam of 263.8 axcmin. Finally, since the real Haslam map 
has been produced combining different surveys and has been 
degraded and reprocessed, it is not straightforward to deter- 
mine the level of pixel noise present in the map, which will 
also be correlated. However, we expect it to be very small 
for the considered scales. Therefore, we have neglected the 
noise when generating the map, which has been sinmlated 
with an effective resolution of 268.7 arcminutes (51 arcmin- 
utes beam coming from the resolution of the Haslam map 
plus a smoothing of FWHM=263.8 arcminutes). 

4.2 Results 

We have applied the method explained in §2 to our simu- 
lated data in order to recover the CMB and galactic fore- 
grounds. Given the resolution and signal to noise of the data, 
we have aimed to reconstruct the different components only 
up to imax = 40, since there is virtually no information 
about the CMB in the data at higher multipoles. We need 
to provide the algorithm with an initial guess for the power 
spectra for each of the components. As shown in H98, the 

* This map is available at http://skymaps.info 
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Figure 4. Input simulated components (from left to right and top to bottom): CMB, free-free, Synchrotron and dust at 50 GHz in 
/iK (thcrmodynamical temperature). The maps have been smoothed with a 7 degree Gaussian beam. The galactic components have 
been plotted in a non-linear scale to allow for a better visualization of the high-galactic latitude region. To allow for a straightforward 
comparison with the reconstructed map, the input components include only information up to £ = 40. This produces the ringing artifacts 
visible in the dust map. 



reconstructions are not very sensitive to this initial guess, 
provided one iterates over the power spectra, i.e. one per- 
forms the reconstruction with an initial guess and uses those 
reconstructed maps to provide starting power spectra for the 
next iteration, until convergence is obtained. For the galac- 
tic components, we have chosen initial power spectra that 
differ appreciably from the original input maps to show the 
performance of the method even when the initial power spec- 
tra are far from the correct ones. Regarding the CMB, we 
have chosen not to iterate over its power spectrum, but to 
start each iteration with a CMB model which is compatible 
(but that differs from the input one) with the power spectra 
derived from the COBE data (at the lowest £'s) and from 
Archeops (at the highest multipoles). As shown in S02, an 
approximated initial guess is enough for MEM to find the 
underlying power spectra. It must also be pointed out that 
if we were using a prior for the CMB power spectrum that 
significantly differs from the true one and we do not iterate 
over this quantity, this could bias the results obtained by the 
method, especially at those scales with a low signal to noise 
ratio. However, we do have a good knowledge of the shape of 
the power spectrum of the CMB at these low ^'s. Therefore 
it is reasonable to make use of these information in order to 
improve the CMB reconstruction. In any case, we would like 
to emphasize that the knowledge of the CMB (or any of the 
other components) power spectrum is not necessary for the 
method to work. If there was a total absence of knowledge 
of the CMB prior, we would need to iterate over its power 
spectrum to find the correct reconstruction. In fact, if we 



were using high quality data, such an analysis without pro- 
viding any prior information should also be performed. This 
would avoid biasing the results as well as point out possible 
inconsistencies with the supplied prior. 

We also need an estimation of the noise for each data 
map. The noise of the COBE maps is well known, but this 
is not the case for the Haslam map, as mentioned in the 
previous section. In order to provide a reasonable value for 
the calculation of the function we have used the following 
trick. The Haslam map has some structure beyond I — 40, 
but we are going to reconstruct the map only up to that mul- 
tipole. Since the map has been smoothed and processed we 
expect the noise to be very small at the scales that we are 
reconstructing. Therefore when subtracting our predicted 
data map (generated using the reconstructions up to imax) 
from the true data map to obtain the x^, the difi'erence will 
come mainly from the power beyond the maximum recon- 
structed multipole and that can be considered our effective 
noise. Therefore, the estimation of the noise for this simu- 
lated data map has been obtained as the dispersion of the 
map obtained subtracting the simulated Haslam map with 
power up to Imax from the same map at full resolution. In 
practice we have used the same trick to estimate the noise of 
the COBE-DIRBE channel. This is due to the fact that the 
noise per pixel of this map is very low after repixelisation 
and smoothing and it was necessary to take into account 
the structure present beyond Imax, which was giving the 
main contribution to the . This estimation of the noise 
produced in the correct range. 
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Table 2. Spectral parameters and value of X^/^/ for the ten best 
cases from low noise simulated data 



Case 


Td 








9 / 

X /"/ 


1 


18 


2.0 


-0.16 


-0.8 


12156 


2 


18 


2.0 


-0.19 


-0.8 


12167 


3 


18 


2.0 


-0.13 


-0.8 


12171 


4 


18 


2.0 


-0.19 


-0.9 


12215 


5 


16 


2.2 


-0.13 


-0.8 


12229 


6 


16 


2.2 


-0.16 


-0.8 


12237 


7 


16 


2.2 


-0.19 


-0.9 


12240 


8 


18 


2.0 


-0.16 


-0.9 


12244 


9 


20 


2.0 


-0.13 


-0.8 


12252 


10 


20 


2.0 


-0.16 


-0.8 


12258 



4-2. 1 Estimation of the spectral parameters in a low noise 
case 

Unfortunately the COBE-DMR data maps arc very noisy, 
so it is difficult to test the performance of our method and, 
in particular, the determination of the spectral parameters, 
with this data set. Therefore, we have first tested our flex- 
ible MEM algorithm in an artificially low noise case. For 
this test, we have used the same simulations described in 
the previous section but the noise level of the COBEl-DMR 
channels has been lowered by a factor of 5. 

We have applied our method to this simulated data 
assuming different sots of spectral parameters. In partic- 
ular, we have used all the possible combinations (a total 
of 81) of the following values: /3ff=-0. 13,-0.16,-0.19, (isyu = 
-0.7, -0.8, -0.9, Td=16,18,20 and ad=l-8,2.0,2.2. For each 
set of spectral parameters we have iterated over the power 
spectra for all the galactic components to find the best pos- 
sible reconstructions in each case. We have then looked at 
the value of the reconstructed maps for each combina^ 
tion of spectral parameters as an indicator of the quality of 
the reconstructions. We find that, as one would expect in 
an ideal case, the reconstruction with the lowest value 
(or equivalently X^/*^/) corresponds to the one obtained us- 
ing the correct set of spectral paramaters (i.e., /3ff=-0.16, 
B«yn = —0.8, Td=18 and ad=2.0). Since the noise of each 
data map is quite small, there is no room for uncertainties 
in the spectral indices and the successfully picks the right 
set of spectral parameters. We have numbered the different 
cases according to the obtained x^ value. Case number 1 cor- 
responds to the lowest x^ whereas case 81 is the one with 
the highest value. Table 2 gives the value of the spectral 
parameters and of x^ /rif for the 10 first cases. 

In addition to the x^ there are other quantities that 
can be calculated from the reconstructed components that 
also provide information about the quality of the recon- 
structions. In particular, we have also studied the behaviour 
of f/uf, S/uci the cross-correlations between the CMB re- 
constructed map and the three galactic components (which 
should be zero) and the dispersion of the CMB reconstructed 
map tJcMB- AH those quantities arc plotted in Fig. 5 versus 
the case number. The error etrue in the CMB reconstruction 
(smoothed with a 7 degree beam) is also shown and has been 
calculated as 

etrue = ^/{{Ti - TrY) - - T.)^ (15) 



Table 3. Dispersion values for input, reconstructed (from case 1) 
and residuals maps smoothed with a 7 degree beam at 50 GHz in 
/iK for the low noise simulated data. These values arc given for 
the whole map (col. 3) , for the region outside the galactic cut 
(col. 4) and for the galactic centre (col. 5). For comparison the 
estimated error is also given in the last column. 



Cpt 


Map 




(Tout 


TGal 


eest 




Input 


36.4 


36.5 


36.3 




CMB 


Rec. 


35.8 


35.4 


36.5 






Resid. 


9.5 


9.4 


9.8 


10.9 




Input 


186.2 


34.8 


255.6 




FF 


Rec. 


186.6 


36.5 


2-56.2 






Kesid. 


11.9 


11. S 


11.9 


2:17 




Input 


42.7 


10.1 


59.2 




Synch. 


Rec. 


42.7 


10.1 


59.2 






Hesid. 


0.55 


0.51 


0.55 


1. 1 




Input 


16.6 


0.66 


25.1 




Dust 


Rec. 


16.6 


0.66 


25.1 






Resid. 


0.08 


0.08 


0.08 


1.3 



where Tj and Tr correspond to the input and reconstructed 
map smoothed with a 7 degree beam. It is interesting to note 
the high correlation between the value of the x^ and the 
CMB reconstruction error. Moreover, this error correlates 
also very strongly with the rest of the plotted quantities. As 
one would expect, the maps with a low etrue also present 
values of the cross-correlations between the CMB and the 
galactic restored maps close to zero. Small values of f> and 
SI correspond as well to small values of etrue. Finally, we 
also find that the lowest values of (Jcmb are also the ones 
with better CMB reconstructions. This can be explained 
taking into account that those CMB recostructions obtained 
with wrong spectral indices will be contaminated by galactic 
emission and therefore the dispersion of the CMB map will 
increase. 

Since the main objective of this test was to study the 
estimation of the spectral parameters using different indi- 
cators of the quality of the reconstructions, we will not go 
into detail with regard to the reconstructed maps and power 
spectra. However we give, as reference, the difference recon- 
structed errors for the four recovered components in Table 3 
for case 1. 

i.2.2 Estimation of the spectral parameters in the 
simulated case with realistic noise 

We have repeated the previous analysis using our (realistic) 
simulated data set and have obtained reconstructions for 
the same 81 combinations of spectral parameters. As in the 
previous case the reconstructions have been obtained iter- 
ating over the power spectra of the galactic components for 
each set of spectral parameters. We have then investigated 
the behaviour of the different quantities that might indicate 
which combination of parameters is correct. If we look at 
the reconstruction with the lowest x? value we find the fol- 
lowing set of parameters: Td = 18, Qd = 2.0, /3ff = —0.13, 
/3syn = —0.7 (corresponding to a X^/^/ = 12191) whereas 
the reconstruction with the correct spectral parameters has 
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Figure 5. The following quantities have been plotted versus the case number for the low- noise sinmlated ease: left column (from top 
to bottom); <p/nf, 'x^ juj, S/ric and dispersion of the reconstructed CMB map. Right column (from top to bottom): cross correlations 
between the reconstructed maps of CMB and free-free, CMB and synchrotron, CMB and dust and error of the CMB reconstruction. For 
reference, the dotted line in the cross correlation panels show the level of the correlation between the input maps of the CMB and the 
corresponding galactic component. 



the second lowest value which is also very close to the 
minimum {y^ jnf = 12192). In fact, 3 more of the recon- 
structions have values ')i jrif < 12193. Therefore we find 
that the data are just too noisy to discriminate clearly be- 
tween different spectral data sets using just the information 
given by the ■ However, we have seen in the previous sec- 
tion that there is more information that we can extract from 
the reconstructed maps to assess the quality of the recon- 
structions and to determine the spectral parameters. There- 
fore, we have also considered the values of the entropy 5", the 
(p function, the dispersion of the reconstructed CMB map 
ffcMB: and the cross-correlation between the reconstructed 
CMB and each of the reconstructed galactic components. In 



particular we have constructed an empirical selection func- 
tion G that is a linear combination of the former quantities 
and is defined as 



G 



ai 1-02 1- as h 04 Cff 

n/ Tie 



+ oslcsynl + aelcdusti + otCcmb 



(16) 



where Cj denotes the cross-correlation between the CMB 
reconstruction and the corresponding galactic component. 
The minimum of G will give us the best reconstructions. 
The chosen weights Oj are given in Table 4 and have been 
determined by looking for an optimal combination that gives 
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Case number Case number 




Case number Case number 



Figure 6. The panels show, for the reahstic simulated data, the residuals dispersion (solid line) versus the case number for each of the 
reconstructed components after smoothing with a 7 degree Gaussian beam. For comparison the estimated error is also plotted (dashed 
line) 



Table 4. Weights for the calculation of G 



Coefficient 


Value 


ai 


0.6 


02 


0.6 


as 


5 


a4 


100 


as 


100 


06 


100 


07 


1 



more weight to those reconstructions with higher quality and 
whose spectral parameters deviate less from the true ones. 
The goal of combining all this information is that even if we 
can not reUably determine the spectral parameters, given 
the low signal to noise of the data, at least we can prevent 
the uncertainties in the knowledge of such parameters from 
severely contaminating the CMB reconstruction. 

Using G, we find that the preferred reconstruction is the 
one with the correct spectral paramaters, whereas that with 
the lowest value is now in position 10 (sec Table 5 to see 
the values of the parameters for the first 10 cases). The dif- 
ferent cases have been ordered according to the value of G, 
so those cases with a lower number correspond to sets of pa- 



Table 5. Spectral parameters and value of G for the ten best 
cases from (realistic) simulated data 



Case 


Ta 




0ff 




G 


1 


18 


2.0 


-0.16 


-0.8 


54.0050 


2 


18 


2.0 


-0.19 


-0.9 


54.0055 


3 


18 


2.0 


-0.13 


-0.8 


54.0063 


4 


18 


2.0 


-0.19 


-0.8 


54.0084 


5 


18 


2.0 


-0.16 


-0.9 


54.0085 


6 


20 


1.8 


-0.19 


-0.7 


54.0525 


7 


18 


2.0 


-0.13 


-0.9 


54.0536 


8 


16 


2.0 


-0.19 


-0.7 


54.0547 


9 


18 


2.0 


-0.16 


-0.7 


54.0556 


10 


18 


2.0 


-0.13 


-0.7 


54.0563 



ramctcrs that produce a smaller G. It is interesting to note 
that the smoothed 7 degree CMB reconstruction error of the 
case prefered by the G quantity (case 1) is 21.5/iK whereas 
that of the case with the minimum (case 10) is 22.9;uK. 
In fact if we calculate the correlation between the values of 
G and of the CMB reconstructed error for the 81 cases we 
find a value of 0.993 versus 0.969 for the correlation between 
the value and the same error. In addition, there is also 
another importan reason to use the combined information 
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Figure 7. The quantities used to calculate G, and G itself, are plotted versus the case number for the simulations. Left column (from top 
to bottom); f/rif, 'x^ juf, S/ric and dispersion of the reconstructed CMB map. Right column (from top to bottom): cross correlations 
between the reconstructed maps of CMB and free-free, CMB and synchrotron, CMB and dust and G. For reference, the dotted line in the 
cross correlation panels show the level of the correlation between the input maps of the CMB and the corresponding galactic component. 



rather than just the x^- For real data, where the foregrounds 
can not be well modelled by a simple law, or where system- 
atics may be present, it seems a good idea to combine all 
the available information. In particular, the will not per- 
form as well as in the ideal case and therefore we should also 
make use of quantities such as the cross-correlations between 
CMB and the galactic components that can be more reliable 
to assess the quality of the reconstruction for real data. 

The behaviour of the quantities used to construct G 
is illustrated in Figs. 6, 7 and 8. The different quantities 
have been plotted versus the case number (ordered accord- 
ing to G); Fig. 6 shows the dispersion of the residuals (solid 
line) for each of the components smoothed with a 7 degrees 
Gaussian beam. Note that in the CMB case, there is a clear 



correlation between the reconstruction error and the order- 
ing of the cases: a lower value of G (which is plotted in the 
bottom right panel of Fig. 7) will produce in general a bet- 
ter CMB reconstruction. It is also striking that this error 
varies very little for the cases with the lowest values of G 
In fact, the difference outside the galactic cut^ between the 
CMB reconstruction (smoothed with a 7 degree Gaussian 
beam) of case 1 and all those up to case 25 is ^ 3/iK, which 
is well below the statistical errors. As expected, the CMB 



t We will always refer to the custom galactic cut of Banday et 
al.(1997) in HEALPix pixelization, which masks a total of 4594 
pixels. 
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Case number 

Figure 8. Spectral parameters versus case number for the simu- 
lations: dust model (top), spectral index for the free-free (mid- 
dle) and spectral index for synchrotron (bottom). The dust 
model number correspond to: T<j=18,ad=2 (1), T^=20,ad=1.8 
(2), Td=16,ad=2 (3), Td=16,ad=2.2 (4), Td=20,ad=2 (5), 
Td=18,ad=2.2 (6), Td=20,ad=2.2 (7), Td=18,ad=1.8 (8), 
Td=16,ad=1.8 (9). 



reconstruction in the galactic region is more dependent on 
the spectral parameters and the differences range between 
3 — 12jj,K for the same cases. The dashed hnes in Fig. 6 in- 
dicate the residuals dispersion estimated by our method. As 
already mentioned, this error is only reliable if the spectral 
parameters arc close to the true values. If this is the case, 
this error gives a good estimation of the residuals dispersion. 

Fig. 7 shows some additional sensitive quantities versus 
the case ordering. We see that low values of x^j V 
the absolute value of the entropy also go in the direction 
of producing better reconstructions. In addition, the cross 
correlations between the reconstructed CMB and galactic 
components are also indicators of the reliability of the re- 
constructions. Another very interesting result is given in the 
bottom left panel of Fig. 7, where the dispersion of the re- 
constructed CMB map smoothed with a 7 degree Gaussian 
beam is given. The correlation between this quantity and 
the actual error of the CMB map is striking. As already 
mentioned, this can be understood since errors in the CMB 
restored map would mainly come from the introduction in 
the reconstruction of galactic contamination, which will give 
rise to a higher dispersion of the map. The fact that different 
foreground models produce similar CMB reconstructions in- 
dicate that we have degenerate cases, due again to the low 
signal-to-noise ratio of our data. 

This effect can be seen in the two bottom panels of 
Fig. 8, where the spectral indices for the free-free and syn- 
chrotron are plotted versus the case number. There is no 
clear trend between those parameters and the quality of the 
reconstruction. Therefore, these data are not precise enough 
to determine unambiguously the free-free and synchrotron 
parameters. However, it is still important to look at all the 



Case number 



Figure 9. Value of G versus case number (ordered from low to 
high value of x^) for the low noise case. 



available information, since not all combinations of 0g and 
Psyn produce equally good CMB reconstructions. In the case 
of the dust component (see bottom right panel of Fig. 6 and 
top panel of Fig. 8) there is a visible correlation between the 
correct dust model and the case ordering. Therefore the dust 
parameters can be determined by the method itself for the 
case of a spatially invariant dust model. In a real data set, 
where spectral variability would occur, the picture would 
not be so clear, but there would still be some trends in the 
graphs that would point out to some preferred models (see 
discussion in next section). 

We note that there are other possible choices of weights 
that would assign the best G to the case with the correct 
spectral parameters. In fact, in practice, it would be very dif- 
ficult to distinguish between the five cases with a lower (and 
almost identical) value of G for this data set, since other sim- 
ilar choices of weights lead to a reordering of the top cases. 
We remark, however, that the CMB reconstruction is very 
robust independently of the chosen spectral parameters for 
those cases with low values of G. The fact that G is so fiat 
indicates again that the data arc too noisy to discriminate 
unambiguosly between the different sets of spectral param- 
eters. MEM is able to accomodate a handful of different 
models within the noise, providing virtually indistinguish- 
able CMB reconstructions. 

It is interesting to look at the result if we apply the 
selection function G to the low noise simulated case of the 
previous section. Fig. 9 shows the values of the G function 
versus the case number in the low noise case. Note that the 
cases are ordered according to the value of the obtained 
in the low noise reconstructions and therefore the numbering 
does not coincide with the one on this section. As expected, 
the G selector works also in this case and the reconstruc- 
tion with the minimum value of G coincides with the one 
obtained with the minimum value of (which had the cor- 
rect spectral parameters). The fact that the value of G is 
higher for the low noise case with respect to the realistic case 
should not be surprising. G can be used to compare different 
reconstructions for the same experiment. However it should 
not be used to compare the quality of reconstructions for ex- 
periments with different characteristics. For instance, MEM 
gives more weight to minimising the value of the x^ than 



© 2003 RAS, MNRAS 000, 000-000 



14 Barreiro et al. 



Table 6. Dispersion values for input, reconstructed (from case 1) 
and residuals maps smoothed with a 7 degree beam at 50 GHz 
in /iK. These values arc given for the whole map (col. 3) , for the 
region outside the galactic cut (col. 4) and for the galactic centre 
(col. 5). For comparison the estimated error is also given in the 
last column. 



Cpt 


Map 




O'out 


<^Gal 






Input 


36.4 


36.5 


36.3 




CMB 


Rec. 


34.3 


33.2 


37.4 






Resid. 


21.5 


21.4 


21.3 


21.8 




Input 


186.2 


34.8 


255.6 




FF 


Rec. 


187.2 


46.9 


256.3 






Resid. 


36.9 


33.2 


42.4 


42.0 




Input 


42.7 


10.1 


59.2 




Synch. 


Rec. 


42.8 


10.2 


59.4 






Resid. 


1.7 


1.5 


1.9 


2.1 




Input 


16.6 


0.66 


25.1 




Dust 


Rec. 


16.6 


0.66 


25.1 






Resid. 


0.04 


0.03 


0.05 


1.3 



that of the entropy in the low noise case producing a higher 
value of S (although with a slightly lower value of a) what 
contributes to a higher G. Also the value of Ucmb is higher 
in the low noise case since the CMB has been reconstructed 
with (intrinsic) higher resolution. A more interesting con- 
clusion can be derived by looking at the range of G for the 
low and realistic noise cases. Whereas G spans from ~ 56 to 
^ 73 for the low noise case, in the realistic case the range 
for the same 81 cases goes from ^ 54 to 56. This shows 
again that the low signal to noise of the realistic data make 
it very difficult to discriminate between the different con- 
sidered cases. We also note that the optimal choice of the 
coefficients defining G may vary depending on the charac- 
teristics of the considered experiment. Therefore, a detailed 
study using simulations should be performed in each case. 

4-2.3 Reconstructed maps and power spectra for the 
simulated case with realistic noise 

We have plotted the reconstructed maps corresponding to 
case 1 (i.e. using the correct spectral parameters) smoothed 
with a 7 degree Gaussian beam in Fig. 10 in the same scale 
as the input maps (Fig. 4) to allow for a straightforward 
comparison. The residuals map for each component is given 
in Fig. 11. The values of the residuals dispersion for each 
component (all-sky, high and low galactic latitude), the dis- 
persion of the input maps and the estimated errors are sum- 
marized in Table 6. 

The dispersion of the CMB residuals on an angular scale 
of 7 degrees is at the level of ~ 21 /iK, which is in good agree- 
ment with the estimated error given by the MEM algorithm 
(see Table 6). Many of the main features of the CMB in- 
put map are also present in the reconstruction. However, 
the smallest structure has clearly been damped in the re- 
constructed map. This is expected since at the highest con- 
sidered i's the COBE-DMR data are dominated by noise 
and MEM just defaults to zero in absence of any useful in- 
formation. It is interesting to point out that the CMB errors 



at high and low galactic latitudes are actually comparable 
(~ 21/iK in both cases, see Table 6), showing that the map is 
equally well recovered independently of the Galaxy. This is 
the case even for small departures of the spectral parameters 
from the true ones, such as those in the 5 best cases, which 
have very similar residuals dispersions in both regions. This 
is due to the very high noise of the COBE-DMR maps that 
allows to accommodate the data to different spectral mod- 
els. In fact, as already mentioned, the CMB reconstruction 
is very robust for models with a low value of G. The differ- 
ence between the 7 degree smoothed CMB reconstruction 
of case 1 and those of cases 2 to 5 are < 2/Ltif outside the 
Galaxy and < AfiK inside the galactic cut for all the cases, 
values which arc significantly lower than the statiscal error. 
Therefore, even if we can not determine with total reliability 
the spectral parameters from our data set due to their low 
signal-to-noise, we can avoid the introduction of errors in the 
CMB reconstruction due to these uncertainties, especially at 
the high galactic latitude region. For cases with larger values 
of G, the reconstructed error inside the Galaxy starts to be 
systematically higher than that of the high galactic latitude 
region. 

Regarding the free-free map, it can be seen that the 
galactic plane is reasonably well recovered (with an error 
~20 per cent) whereas most of the structure outside the 
galactic cut has been lost. In fact much of the signal re- 
covered at high galactic latitude takes negative values. The 
synchrotron and dust maps have been very well recovered 
since MEM has succeeded in tracing these emissions from 
the Haslam and COBE!-DIRBE maps respectively. There are 
only some small differences between the input and recon- 
structed maps, mostly at small scales, and the dispersions 
of the residuals are at the level of ~ 4 per cent for the syn- 
chrotron and ~ 0.3 per cent for the dust component. 

Of course, when the wrong spectral dependence is as- 
sumed, the errors increase appreciably for all the galactic 
components. However, these differences come mainly from a 
normalization factor rather than from the spatial structure. 
For instance the reconstruction error for the synchrotron in 
case 2 (where /3syn = —0.9 was assumed) is ^ 17/^K (ver- 
sus ~ 2/iK for case 1) whereas the spatial cross correlation 
between the input and reconstructed smoothed maps is at 
the same level (0.999) for both cases. Thus, MEM is find- 
ing the right amplitude and structure of the synchrotron at 
the Haslam frequency, where this emission dominates, and 
then extrapolating to the reference frequency using the con- 
sidered spectral index. Similar idcEis apply to the dust and 
free-free emissions. 

In Fig. 12 the true (thin solid line) and reconstructed 
(thick solid line) power spectra for the unsmoothcd input 
and recovered maps have been plotted. The dotted lines are 
the estimated la confidence level for the reconstructed com- 
ponents. Finally, the dashed line correspond to the initial 
power spectra supplied to the MEM algorithm. In the CMB 
panel (top right) we have also plotted the power spectra 
measured from COBE (solid squares). Note that the CMB 
power spectrum used as initial guess in MEM (dashed line) 
is similar but differs from the true one (thin line), allow- 
ing for possible errors in the estimation of the prior. Even 
without iterating in the CMB power spectra, the recovered 
Ci's follow quite well the true ones up to £ ~ 15 and then 
start to drop due to the resolution of the COBEl-DMR data. 
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Figure 10. Reconstructed components from realistic simulated data at 50 GHz obtained using the correct spectral parameters (case 1). 
The maps have been smoothed with a 7 degree Gaussian beam and have been plotted in the same scale as the input maps of Fig. 4 to 
allow a straightforward comparison 
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Figure 11. Residuals of the simulated data for the correct spectral parameters (case 1) obtained by subtracting the smoothed input 
maps (Fig. 4) from the smoothed reconstructions (Fig. 10). All the maps have been plotted in linear scale. 
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Figure 12. Original (thin solid), initial guess (dashed) and reconstructed (thick solid) power spectra for the (realistic) simulated data. 
The dotted lines indicate the estimated Itr confidence level for the reconstructed power spectra. The correct spectral parameters have 
been used for the reconstructions (case 1). For comparison the power spectra measurements obtained from the COBE-DMR data are 
also shown (solid squares in top left panel, taken from Tegmark & Hamilton 1997). 



As already mentioned, this reflects a loss of resolution in 
the reconstructed map with respect to the input one. We 
may also wonder about the quality of the recovered power 
spectra inside and outside the galactic cut. Fig. 13 shows 
the true and reconstructed power spectra in these two re- 
gions of the sky. For simplicity, we have estimated the Ce's 
padding with zeros the masked region and they have been 
rescaled according to the considered area. For the CMB case 
the results are quite similar for both regions of the sky, and 
the reconstructed power spectrum follows approximately the 
true power up to f ~ 15, as it was the case for the whole 
sky. 

Regarding the galactic components, it is quite striking 
that MEM is able to recover the right power spectra even 
when the initial guess was far off from the correct Ce's. The 
free-free power spectra (top right panel of Fig. 12) has been 
reasonably well recovered up to ^ ~ 30 (although some ex- 
cess is present at the smallest scales) where it sharply drops 
to zero. As in the case of the CMB, the free-free informa- 
tion comes mainly from the noisy COBE-DMR data whose 
signal to noise ratio and resolution does not allow MEM to 
provide a reconstruction at higher i's. An interesting result 
is also found by looking at Fig. 13. The reconstructed FF 
power spectrum inside the galactic plane follows quite well 
the input one up to / 25, however the reconstruction is 



much poorer outside the galactic cut. In fact, a clear excess 
of power is seen at ^ ^ 10, which gives rise to the spuri- 
ous structure that is present in the FF reconstructed map 
at high galactic latitude and to the excess of power seen in 
the whole-sky reconstructed power spectrum. The reason for 
MEM to fail in recovering the FF in this region of the sky is 
again the low signal-to-noise ratio of the data together with 
the weakness of the FF component outside the Galaxy. The 
synchrotron is faithfully recovered up to ^ ^ 20, result that 
holds for both the high and low galactic latitude regions as 
well as for the whole map reconstruction. At higher ^'s there 
is some excess power which is responsible for the structure 
seen in the synchrotron residuals. Finally, the recovered dust 
power spectra in the galactic region follows very well the in- 
put one up to practically the considered imax- This is due 
to the fact that the dust emission is mainly recovered from 
the COBE-DIRBE map, which provides the algorithm with 
enough resolution to recover faithfully this component up 
to very high £. Regarding the high galactic latitude region, 
there is an excess of power present at ^ ^ 20. However this 
structure does not show in the global power spectrum re- 
construction, which is very good up to high £, since this is 
dominated by the emission at the galactic region, which has 
been very well recovered. 
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Figure 13. Original (thin solid) and reconstructed (thick solid) power spectra for the simulated data inside and outside the galactic 
plane. The correct spectral parameters have been used for the reconstructions (case 1). 



5 ANALYSIS OF REAL DATA 

We have applied our MEM algorithm to the set of data de- 
scribed in §3. As in the previous section, we have aimed to 
reconstruct the CMB, free-free, synchrotron and dust emis- 
sions using the COBE-DMR, Haslam and the lowest fre- 
quency COBE-DIRBE maps. As the initial guess for our 
reconstructions we have chosen the power spectra of the sim- 
ulated components of the previous sections (smoothed with 
the COBE-DMR beam for the CMB and the free-freet and 



t The two all-sky templates for free-free emission (Dickinson et 
al. 2003, Finkbeiner 2003) became available shortly after the com- 
pletion of this work. Therefore, instead of our partially mock tem- 
plate, we could have used one of this new maps to compute the 
initial free-free power spectrum. However, as shown using simu- 



with a Gaussian beam of FWHM = 263.8 arcmin for the 
dust and synchrotron). As before, we have iterated over the 
power spectra for the galactic components but not for the 
CMB, since we have some prior knowledge about its power 
spectra that we would like to include in the analysis. 



5.1 Estimation of the spectral parameters 

First of all we have applied MEM for 180 different sets of 
spectral parameters with values in the range 16 to 22 K 



lated data, the final result is very insensitive to the initial guess 
for the power spectrum. Thus we do not expect that our results 
would be appreciably modified by using a more realistic free- free 
template. 
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Table 7. Spectral parameters and value of G for the ten best 

cases from real data 
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Case number 
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Figure 14. Dust model (top), free- free spectral index (mid- 
dle) and synchrotron spectral index (bottom) versus the case 
number for real data. The dust model number correspond to: 
Td = 19,ad = 2 (1), Td = 18,ad = 2 (2), = 16,ad = 2.2 
{3),Ta = 22,ad = 1.8 (4), = 20,ad = 2 {5),Ta = 16,ad = 2 
(6),Td = 20,ad = 1.8 (7), = 19,ad = 1.8 (8),Td = 22,ad = 2 
(9),Td = 18,ad = 2.2 (10), Ta = I8,aa = 1.8 (ll),Td = 
20,ad = 2.2 (12),Td = 19,ad = 2.2 (13),Td = 22,ad = 1.6 (14), 
Td = 16,ad = 1.8 (15),rd = 20,ad = 1.6 (16), Tj. = 18,ad = 1.6 
(17),Td = 16,ad = 1.6 (18) 



Case 


Td 




l^ff 


Psyn 


G 


1 


19 


2.0 


-0.19 


-0.8 


54.6674 


2 


18 


2.0 


-0.19 


-0.8 


54.6678 


3 


18 


2.0 


-0.19 


-0.7 


54.6747 


4 


16 


2.2 


-0.19 


-0.9 


54.6812 


5 


22 


1.8 


-0.19 


-0.7 


54.6933 


6 


20 


2.0 


-0.19 


-0.9 


54.6946 


7 


18 


2.0 


-0.19 


-0.9 


54.7056 


8 


19 


2.0 


-0.16 


-0.8 


54.7216 


9 


18 


2.0 


-0.16 


-0.7 


54.7223 


10 


22 


1.8 


-0.19 


-0.8 


54.7227 



Table 8. Dispersion values for the reconstructed components 
from real data smoothed with a 7 degree Gaussian beam at 50 
GHz in //K (for case 1) for the realistic simulated data. These 
values axe given for the whole map (col. 2), for the region outside 
the galactic cut (col. 3) and for the galactic centre (col.4). The 
estimated error is also given in column 5. 



Component 




(Tout 




Best 


CMB 


41.1 


35.0 


47.4 


21.8 


FF 


238.4 


40.9 


360.7 


47.7 


Synch. 


33.6 


10.7 


45.6 


2.6 


Dust 


12.6 


0.67 


18.6 


0.88 



for the dust temperature, 1.6 to 2.2 for the dust emissivity, 
-0.9 to -0.6 for I3,yn and -0.19 to -0.13 for /3//. We have 
then ordered the different cases according to the obtained 
value of G. The lowest value of G is obtained for the param- 
eters Td = 19 K, Ud = 2, Pff = -0.19 and Psyn = -0.8. 
Fig. 14 shows the values of the spectral parameters versus 
the case number. In Fig. 15 the value of G as well as the 
quantities used to calculate it are given for the 180 cases. 
It can be seen that the trends of the curves are the same 
as those of the simulated cases (compare with Fig. 7). Not 
surprisingly, the value of (p, x^-i \S\ and ctc'mb are now a 
bit higher than in the simulated case, which indicates that 
our model does not perfectly fit the data. This is due to 
the large number of uncertainties present in the real case, 
such as the spectral dependence of each component, its po- 
sition and/or frequency variability or even the number of 
components. However, the is still very reasonable since 
it is of the order of the number of data. As in the case of 
the simulated data, there is a clear correlation betwoou the 
dust model (top panel) and the value of G, indicating that 
some of the studied dust parameters are preferred by the 
data. In particular models with Od = 2 and Td around 18- 
20 K are in the first positions. Conversely, the considered 
models with = 1.6 are clearly disfavoured. As expected, 
the distinction between the models is less clear than for the 
simulated data, since degenerate cases are more important 
for true data than for ideal ones. Regarding the free-free 
spectral index (middle panel), there is a slight trend from 
the data to favour those cases with more negative values 
of /?//. Values of /3syn between —0.9 and —0.7 are basically 



equally preferred by the data, whereas l3sy„ — —0.6 seems 
to produce reconstructions with higher values of G. Finally, 
for comparison with the ideal simulated case, we have plot- 
ted the estimated error for the reconstructed components in 
Fig. 16. They are at similar levels as those in the simulated 
cases. 

Table 7 gives the ten cases with the lowest values of G. 
The difference between the CMB reconstruction for the se- 
lected case one and all the cases up to number 15 is between 
0.8 and 2.8 /uK outside the galactic cut and between 1.6 
and 10 in the galactic centre. Therefore, as happened in 
the ideal case, the CMB recovered map in the high galactic 
latitude region is very robust against certain variations of 
spectral parameters. 

5.2 Reconstructed maps and power spectra 

Fig. 17 shows the reconstructed CMB, free-free. Synchrotron 
and dust cnussions for case 1 {Td = 19, Ud = 2, Pff = —0.19 
and (3syn = —0.8). 

Table 8 gives the dispersion level of the reconstructed 
components smoothed with a 7 degree Gaussian beam at the 
reference frequency of 50 GHz, for the whole-sky, the Galaxy 
and the high galactic latitude region. For the CMB, we can 
see that the dispersion for the two considered regions differ 
appreciably (35 versus 47 /iK). This indicates that there is 
some contamination at low galactic latitudes in the CMB 
reconstructed map. This is not surprising given the uncer- 
tainties present in our model of the data (spatial variability 
of spectral parameters, frequency model for each component. 
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Figure 15. As Fig. 7 but for the 180 cases used to recover the microwave sky from the real data. 



number of present components, possible systematics, etc.). 
The estimated error of the CMB reconstruction on a scale of 
7 degree is ~ 22/iK. However, given the contamination of the 
CMB map in the galactic centre, we can consider that this 
value is a fair estimate of the CMB residuals only outside 
the galactic cut. 

For comparison we have also plotted in Fig. 18 a CMB 
map obtained by coadding the 53 and 90 GHz COBE-DMR 
maps (each pixel weighted according to the inverse of its 
noise variance) which has been smoothed with a 7 degree 
Gaussian beam. To allow for a better comparison, we have 
removed the monopole and dipole outside the galactic cut 
which were mainly due to galactic emission. As expected, we 
find this map to be well correlated with our reconstructed 
CMB map (whose monopole and dipole have also been re- 
moved outside the galactic mask to calculate this correlar 
tion), at the level of 0.85 (outside the galactic cut), which 



confirms the presence of CMB signal in our reconstruction. 
The difference between this coadded map and our recon- 
struction outside the galactic cut is also shown in the same 
figure (bottom). Note the smaller scale range in the differ- 
ence map which indicates the presence of common structure 
in both maps. However, some clear structure is also seen 
in the difference map. This is not surprising since not at- 
tempt to remove foreground emission has been made in the 
coadded map. 

Regarding the free-free reconstructed map (top right 
panel of Fig. 17), we are able to recover emission only inside 
the galactic cut. As in the simulated case, the free-free signal 
at high galactic latitude is lost. In fact, the signal recovered 
outside the galactic cut takes mainly negative values and 
its dispersion is lower than the estimated error (~ 48 ijK). 
Regarding the galactic centre, MEM has only been able to 
recover emission in ~ 30 per cent of the area of the galactic 
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Figure 16. Estimated errors of eacii of tfie reconstructed components versus tiie case number for tfie real data. Each case corresponds 
to a diflferent set of spectral parameters and they have been ordered according to the value of G. 



cut, whereas the rest of the signal has zero or negative val- 
ues. The dispersion level of the 7 degree smoothed free-free 
component in the galactic plane is ~ 238/iK considering all 
pixels and it increases to ~ 491/iK if we consider only the 
fraction of pixels with physical (i.e. positive) emission. This 
is several times higher than expected from current Ha esti- 
mations and raises again the issue of an anomalous compo- 
nent. It is also interesting to note the high spatial correlation 
between the galactic plane of the reconstructed free- free map 
and that of the COBE-DMR frequency channels. 

The synchrotron reconstructed map (bottom left panel 
of Fig. 17) presents structure both inside and outside the 
Galaxy. As expected, the reconstructed emission is very well 
correlated with the Haslam map (0.96 for the smoothed 
maps). We find a level of ~ 46/nK for the synchrotron emis- 
sion in the Galaxy and of ~ 11/iK outside the Galaxy. The 
estimated error in the reconstruction is 2.6/iK. We should 
note, however, that this is the statistical error associated to 
our method but it docs not take into account the uncertain- 
ties in the determination of /9syn. 

Finally, the bottom right panel of Fig. 17 shows the re- 
constructed dust emission smoothed with a 7 degree Gaus- 
sian beam. The visible ringing is an artifact due to the fact 
that wc arc recovering the dust only up to £max = 40, 
whereas the DIRBE data map (from which the algorithm 
traces this component) has power at higher multipoles. We 
have found a dispersion value for the dust component of 



~ 19^K and of ~ 0.7/iK inside and outside the Galaxy re- 
spectively at 50 GHz. The (statistical) estimated error is 
below 1 nK. 

The (unsmoothed) reconstructed power spectra (solid 
line) for the four recovered components are given in Fig. 19. 
The dashed lines show our initial guess for the power spectra 
and the dotted lines indicate the la confidence level for the 
reconstructed power spectra. Fig. 20 shows the power spec- 
tra for the four reconstructed components inside and outside 
the galactic cut, that have been estimated in the same way 
as for the simulations. 

For comparison we have also plotted the CMB measure- 
ments obtained for COBE-DMR in the CMB panel (top-left) 
of Fig. 19. It is interesting to note that the recovered C'e's 
follow the shape of the power spectra measured from the 
COBE-DMR data but with a higher normalization. In ad- 
dition, the reconstructed CMB power spectra for the high 
galactic latitude region (top left panel of Fig. 20) has ap- 
proximately the expected amplitude whereas the one at the 
Galaxy (top right panel of Fig. 20) presents an excess of 
power at all scales. These arc also indications of the fact that 
we have some galactic contamination in our reconstructed 
map, especially at high galactic latitudes. 

The free-free reconstructed power spectra presents the 
same behaviour as in the simulated case. It seems to be 
recovered up to £ = 30 and then drops sharply to zero, due to 
the low resolution of the COBE-DMR channels, which MEM 
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Figure 17. Recovered components smoothed with a 7 degree Gaussian beam from real data for the spectral parameters T(j=19 K, o = 2, 
Pff = —0.19 and Psyn = —0.8 (case 1) at 50 GHz in fiK. The synchrotron and dust components have been plotted in a non-linear scale. 



has mainly used to recover the free-free signal. We also find 
in the reconstructed power spectra inside and outside the 
Galaxy the same pattern as in the simulations. In particular, 
it seems to be an excess of power at large I at high galactic 
latitude, which produces all the spurious signal that is found 
in this region of the FF reconstructed map. 

The shape of the synchrotron reconstructed power spec- 
tra (bottom left panel) follows quite well that of the input 
one up to ^ ~ 25 but with a lower normalization. This is 
expected since the initial power spectra has been obtained 
from a smoothed extrapolation of the Haslam map up to 
the reference frequency using a spectral index of —0.8 (the 
same value as /Sayn in our case). This extrapolation provides 
therefore an upper limit for the reconstructed synchrotron 
emission (since the Haslam map is dominated by this compo- 
nent but also contains contributions from other emissions), 
which is consistent with our result. At higher ^'s, the power 
spectrum starts to oscillate wildly due to the lack of infor- 
mation in the data to recover this emission. 

The bottom right panel shows the reconstructed Ci's 
for the dust component. As in the case of the simulations, 
the power spectra seems to be well recovered up to ~ imax- 
As expected it follows the shape of the initial guess (the 
differences at the higher ^'s are due to the fact that the input 
has been convolved with a Gaussian beam of FWHM=263.8 
arcminutes) which has been obtained by interpolating the 
COBE-DIRBE channel down to 50 GHz using a dust model 
with Td = 18 and Ud = 2 (different from the best model 
found for the reconstructions). As in the simulated case, 
we also find what seems to be an excess of power at high 



£ values for the reconstructed power spectrum outside the 
galactic plane. 



6 DISCUSSION AND CONCLUSIONS 

We have presented a flexible MEM algorithm that combines 
the advantages of real (or data) space and harmonic space. 
On the one hand, the i^ calculated in data space, which 
allows one to include straightforwardly the properties of the 
noise as well as incomplete sky coverage. On the other hand, 
the entropy is estimated in harmonic space allowing easy 
introduction of available prior information about the power 
spectra of the components that we aim to reconstruct. In ad- 
dition, the method takes into account correlations between 
different modes because we perform a single minimisation 
of ~ nc£^ax variables instead of minimising mode-by-mode 
as in harmonic MEM. We perform this global minimisation 
using a Newton-Raphson method, which needs the Hessian 
matrix H to be evaluated and inverted. In practice, this is 
not feasible since we have a very large number of variables so 
we have approximated H by a block diagonal matrix, which 
gives very good results. Unfortunately, this minimisation is 
very time consuming and the method is many times slower 
than harmonic MEM. 

To test the performance of the method we have ap- 
plied it to simulated spherical data and then to real data. 
In particular, we have used the three frequency channels 
of COBE-DMR, the Haslam map and the lowest frequency 
map of COBE-DIRBE to reconstruct the CMB, free-free, 
synchrotron and dust emissions. 
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Figure 19. Initial guess (dasiied line) and reconstructed (solid line) power spectra for each of reconstructed components obtained from 
the real data for case 1 at 50 GHz. The dotted lines correspond to the Itr confidence level of the reconstructed power spectra. For 
comparison the power spectrum measurements obtained from the COBE-DMR data are also plotted in the CMB panel. 



6.1 Analysis of simulated data 

An important issue that has been thoroughly studied in the 
present work is the determination of the frequency depen- 
dence of the different components. We have modelled the 
free-free and synchrotron with a power law parametrised 
by indices /3ff and /3ayn respectively and the dust by a grey 
body model with two parameters (the dust temperature 
Td and the dust emissivity Od)- We have then applied the 
MEM algorithm to our simulated data (generated using 
/3ff = -0.16, /3syn = -0.8, Td = 18 and ad = 2) for many 
combinations of the spectral parameters and have studied 
the behaviour of different diagnostic quantities when the as- 
sumed values of these parameters depart from the true ones. 

Firstly, we have studied a low noise case where the noise 
level of the COBE-DMR simulated maps have been lowered 
by a factor of 5. In this case we find that the value of 
the reconstructions can successfully identify the correct set 
of spectral parameters. In addition, we also find that there 
exists a high correlation between the error of the CMB re- 
construction and other quantities that can be directly esti- 
mated from the reconstructed maps. In particular, together 
with the value, the entropy, the if function , the dispersion 
of the CMB reconstructed map, and the cross correlations 
between the reconstructed CMB and galactic components 
are good indicators of the quality of the reconstructions. 



Secondly, we have applied MEM to our set of simulated 
data using realistic levels of noise for the same sets of spec- 
tral parameters. Due to the low signal to noise ratio of our 
data set, we find that the is not longer able to pick the 
right combination of spectral parameters. Therefore we have 
constructed an empirical selection function G which is given 
by a linear combination of x^, \S\, ^ , the dispersion of the 
CMB reconstructed map, and the cross correlations between 
the reconstructed CMB and galactic components. We pick as 
our best set of reconstructions the one with the lowest value 
of G. We have considerd a total of 81 different cases, with 
three possible values for each spectral parameter. Using the 
quantity G, we can clearly see that the correct dust model is 
clearly preferred over the others (see Fig. 8). However, it is 
very difficult to pick the correct values of the spectral indices 
for the free- free and synchrotron, since our data do not have 
enough information. In fact the quantity G is almost iden- 
tical for the 5 best cases due to the fact that our noisy data 
can accomodate a range of different spectral models. In any 
case, even if we can not determine unambiguosly the spec- 
tral parameters for the free-free and synchrotron component, 
the CMB reconstruction is very robust outside the galactic 
cut and the differences between the CMB reconstructions 
for the cases with lower values of G are well within the sta- 
tistical errors. In particular the level of the dispersion of the 
residuals of the smoothed 7 degree reconstructed CMB map 
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Figure 20. Reconstructed power spectra for case 1 inside and outside tiie galactic cut for real data. 



outside the galactic cut ranges between 21.2 and 22.0 /iK for 
all the cases between 1 and 30. To illustrate this point we can 
also look at the differences in the CMB reconstructed power 
spectra for the different cases. Fig. 21 shows the CMB power 
spectrum obtained for case 1 (solid line) versus the average 
power spectrum obtained from the 5 cases with lower val- 
ues of G. The error bars correspond to the dispersion in the 
values of the Ce's obtained from the former 5 cases. Note 
that this dispersion is very small, which shows that the re- 
covered CMB power spectrum is also very robust for these 
5 cases. This also indicates that, at least for ideal simulated 
foregrouns with spatially constant spectral parameters, the 
error introduced in the CMB power spectrum due to the 
wrong identification of these parameters (provided we pick 
one of the reconstructions with a lower value of G) is much 
smaller than the statistical error. 

For case 1 (the case with the correct spectral parame- 



ters) we have found that the CMB error reconstruction for 
the (realistic) simulated data is at the level of 21 ^K. The 
level of CMB residuals are comparable inside and outside 
the galactic cut, which indicates that the CMB has been 
equally well recovered in both regions of the sky. Note that 
the estimated error given by the method is ~ 22fiK. The 
recovered CMB power spectrum follows the input one up 
to ^ ~ 15 and then starts to drop due to the resolution of 
the COBE-DMR data. This result is again comparable in- 
side and outside the galactic plane. The free-free emission 
has been recovered with a residual error of ~ 42fiK in the 
galactic plane at 50 GHz. At high galactic latitude, the free- 
free reconstruction has basicly been lost. This result also 
shows in the reconstructed power spectrum, where we find 
a clear excess of structure at ^ ^ 10 outside the galactic cut, 
whereas MEM is able to recover it up to ^ ~ 25 inside the 
Galaxy. The synchrotron and dust emissions have been very 
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imax for the dust. However, the recovered dust power spec- 
trum shows some excess at small scales when only the region 
outside the galactic cut is considered. Note that MEM has 
been able to recover the power spectra of the galactic com- 
ponents, independently of the supplied initial power spectra 
(dashed lines in Fig. 12), which was very far from the true 
one. We have tested that different initial power spectra lead 
to very similar reconstructions, provided MEM uses the it- 
erative mode to obtain the correct power spectra. 




Figure 18. Top: coadded 53 and 90 GHz COBE-DMR maps 
smoothed with a 7 degree Gaussian beam in fjK outside the galac- 
tic cut. The map is plotted in the same scale as the CMB recon- 
struction of Fig. 17. Bottom: the difference between the top map 
minus the CMB reconstructed map of Fig. 17 outside the galactic 
cut is shown. For the comparison, the monopole and dipole out- 
side the galactic region have been subtracted from the coadded 
and reconstructed maps. 
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Figure 21. CMB reconstructed power spectrum for case 1 of the 
simulated data (solid line) and average CMB reconstructed power 
spectra for the 5 cases with a lower value of G (dashed line) . The 
error bars are the dispersion obtained from these 5 power spectra. 

well recovered in the galactic plane with errors of ~ 3 per 
cent and ~ 0.2 per cent, respectively. Outside the galactic 
cut the reconstruction is also quite good, at the level of ~ 15 
per cent for the synchrotron and ~ 5 per cent for the dust. 
The whole-sky power spectra follows also quite faithfully 
the input one up to ^ ^ 30 for the synchrotron and up to 



6.2 Analysis of real data 

In §5 we have applied the same method to the true data. We 
have calculated the value of G for a total of 180 different sets 
of spectral parameters in the range —0.19 ^ /3ft ^ —0.13, 
-0.9 ^ /3syn 5? -0.6, 16 s; Td s; 22 and 1.6 < ^ 2.2. The 
best reconstructions are found for l3g = —0.19, /3syn = —0.8, 
Td = 19 and ad ~ 2. The cross correlations between the 
CMB and galactic components are at a similar level to those 
of the simulations (Fig. 15 versus Fig. 7). However, the value 
of (f, x^, \S\, o"cMB Bind thus of G are a bit higher than the 
ones found for the simulated cases. This indicates that our 
model does not fit the data as well as in the ideal case. This 
is not surprising taking into account the large number of 
uncertainties of the data: the frequency dependence of the 
model, variation of the spectral parameters across the sky 
or in the frequency range or even the presence of some un- 
known component. In spite of this, we still have a reasonable 
fit to our data (the reduced is ~ 1). As in the simulated 
case, some dust models are clearly favoured with respect to 
the others (Fig. 14) although the situation is not so clear as 
in the ideal case. Again this is due to all the uncertainties 
we have with respect to the frequency behaviour of the com- 
ponents which make the presence of degenerate cases even 
more important than in the simulated case. 

The reconstructed maps for the four components are 
given in Fig. 17 for case 1 (the one with the lowest values of 
G) and the corresponding dispersions of the reconstructions 
as well as the estimated errors are shown in Table 8. From 
these results, we can see that there is some clear contami- 
nation of the CMB reconstruction in the galactic plane. For 
instance the dispersion of the recovered CMB is at the level 
of 35 outside the galactic cut, whereas in the galactic 
centre has a value of ~ 47/xK. This is another indication of 
the fact that our model for the galactic components is not 
properly fitting the data. However, we believe that the re- 
constructed CMB sky outside the Galaxy is reasonably well 
recovered. Even if our assumptions about the spectral be- 
haviour of the components are wrong, we have seen in the 
simulations that this had very little effect in the CMB re- 
constructed map outside the galactic cut, providing we are 
taking cases with low values of G. This is also happening 
when using real data. In particular, the difference between 
the 7 degree smoothed CMB restored maps at high galactic 
latitude for case 1 and all cases up to case 20 is ^ 3/xK. 
For the same cases, the differences are ^ 10/iK in the galac- 
tic plane, except for one of the cases. We have also plotted 
the average and the dispersion of the CMB recovered power 
spectrum for the 5 cases with a lower value of G (Fig. 22). 
As in the simulated case, we find very little dispersion within 
these reconstructed power spectra, which confirms that the 
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Figure 22. CMB reconstructed power spectrum for case 1 of 
the real data (solid line) and average CMB reconstructed power 
spectra for the 5 cases with lower values of G (dashed line). The 
error bars are the dispersion obtained from these 5 power spectra. 

CMB reconstruction is not affected by certain variations on 
the spectral parameters. 

As in the simulated case, the free-free emission has 
been recovered only in the galactic centre with a disper- 
sion value at 50 GHz of ~ 238/iK and with an estimated 
error of ~ 48/iK. This value is several times higher than ex- 
pected from Ha measurements and raise again the issue of 
an unknown component. The synchrotron and dust recon- 
structions, however, look very good with estimated errors 
of ~ 3 and ~ 1/xK respectively. These values account only 
for the statistical error but do not take into account the 
uncertainties coming from the determination of the spectral 
parameters. 

As shown in Fig. 19 the CMB recovered power spectrum 
follows the shape of the measurements obtained from the 
COBE data, but with a slightly higher normalisation which 
indicates again the presence of some excess contamination in 
the recovered CMB map. This is confirmed by looking at the 
reconstructed CMB power spectrum inside and outside the 
galactic cut (Fig. 20). At the high galactic latitude region 
the amplitude of the C'l's is consistent with the expected 
level whereas those obtained for the galactic plane present 
an excess of power. The free-free is recovered up to ^ ^ 30 
and then drops to zero, since the data do not have enough 
information at those scales. As in the simulated case, what 
seems a spurious excess of power is found at high ^'s in the 
reconstructed power spectrum of the high galactic latitude 
region. 

The recovered synchrotron power spectrum follows ap- 
proximately the shape of the input power spectra (as ex- 
pected since it has been obtained extrapolating the Haslam 
map) up to ^ 25 . For higher multipoles, the power spectra 
starts to oscillate which indicates that the reconstruction is 
not reliable at those scales. Finally, the dust emission seems 
to be well reconstructed up to the considered imax for the 
whole sky power spectrum. However, some excess is found 
at high £'s in the power spectrum obtained from the high 
galactic latitude region. We would like to point out that, 
as in the simulated case, changing the initial guess power 
spectra has very little effect in the reconstructions. 




Figure 23. Normalised residual noise for the COBE-DMR 53 
GHz channel for case 1 of simulated data (top) and of real data 
(bottom). 

An extra check of the quality of the reconstructions can 
be performed by looking at the normalised residual noise 
for each frequency channel. This map is constructed by sub- 
tracting the predicted noiseless data (obtained using the re- 
constructed maps) from the actual noisy data. This gives 
the 'predicted noise' at the given frequency channel, which 
is then divided by the noise dispersion at each pixel. If the fit 
is acceptable this map should be a realisation of a Gaussian 
white noise field of unit dispersion with no visible struc- 
ture, although one should be careful since incorrect choices 
of the spectral parameters can also give rise to this result 
due to degeneracies. As an illustration, the top panel of 
Fig. 23 shows the normalised residual noise for the 53 GHz 
frequency channel for case 1 of the simulations. No visible 
structure is present in the map, which confirms on the one 
hand, that the anisotropic noise has been properly taken 
into account and, on the other hand, that the Galaxy has 
been well fitted by our components. The bottom panel of the 
figure shows the same map for case 1 of the real data, which 
shows some visible structure in the galactic centre. This fur- 
ther confirms the fact that we are not fully subtracting the 
Galaxy. 

6.3 Two-components dust model 

This inability of the components to fit the data may be due 
either to variations of the spectral parameters with position 
or the frequency range, or to the fact that some other un- 
known component is needed to fit properly the data. A point 
raised by Jones et al. (1999) is the possiblity of reconstruct- 
ing more than one dust component to account for spatially 
variations of the dust spectral parameters. Although each 
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Figure 24. Power spectrum obtained from the (unsmoothed) 
CMB reconstruction of Fig. 17 (solid line) versus the power spec- 
tra of the recovered CMB obtained using two dust components 
(dotted line). The power spectrum measurements obtained from 
the COBE-DMR data are also plotted (solid squares). 

of the dust reconstructions may not have physical meaning, 
adding them together improved the reconstruction. More- 
over, when a dust component with emissivity variations from 
pixel to pixel was included in the simulated data, the er- 
rors of the CMB reconstruction were also reduced by using 
three dust components with different emissivities, as com- 
pared with the case of reconstructing a single dust compo- 
nent whose emissivity was the average of the input template. 

We have tested this possibility by adding to our set of 
data the second COBE-DIRBE channel (at 2141 GHz) and 
reconstructing two dust components. Since we have a total 
of 6 spectral parameters, it is not possible to cover exhaus- 
tively the full range of possible values. Instead, we have fixed 
the values of Ph, /3syn and Od to the ones found to be opti- 
mal in §5.1 (—0.19, —0.8 and 2 respectively) and varied the 
dust temperatures of the two dust components in the range 
of 4 to 22 K. The best reconstruction was found for Ti = 20 
K and r2 = 11 K with a value of G = 53.9005. The CMB 
recovered map is given in the top panel of Fig. 25 and has 
been plotted in the same scale as the one of Fig. 17 to allow 
for a straightforward comparison. The middle panel shows 
the difference between the reconstructed CMB map using 
one dust component minus the one obtained with two dust 
components. It is obvious that the dispersion of the CMB 
map has been lowered inside the Galaxy (47.4 /iK in the one 
dust component versus 42.4 in the two dust components) 
and the new CMB reconstruction seems to have less galac- 
tic contamination. This can also be seen in Fig. 24 which 
shows the power spectra of both reconstructions as well as 
the COBE-DMR data points. The CMB reconstruction ob- 
tained using two dust components is closer to the measured 
power spectra. However, both maps are very similar outside 
the galactic cut. The dispersion of the difTerence map in this 
region is 5 ^K, which is much lower than the estimated sta- 
tistical error. The free-free and synchrotron reconstructed 
maps are very similar for both models, with differences ~3 
per cent for the free- free and ~2 per cent for the synchrotron. 
The dust emission is now composed of the sum of two dust 
components. Each of them may not have a physical meaning 



but the extra degree of freedom should help to account for 
temperature variations of the dust component. In fact, the 
cold dust component presents negative as well as positive 
values in the reconstructed map and should be understood 
as a correction to the hot component, which helps to fit 
both of the COBE-DIRBE maps at high frequencies as well 
as some excess emission at lower frequencies. The total dust 
emission at 50 GHz is given in the bottom panel of Fig. 25. 
Note that the emission in the Galaxy centre has been ex- 
tended with respect to the case when only a single dust 
component was used. Basically, the excess of emission of the 
CMB reconstructed map in the one single component case 
is being put into this cold dust component. Although this 
seems to improve the results, the normalised residual noise 
maps for the COBE-DMR frequency channels present some 
structure in the galactic centre and they look very similar 
to those obtained for the case of one single dust component 
(bottom panel of Fig. 23) . Therefore, our model is still not 
fitting the Galaxy properly and a more exhaustive analysis 
should be performed in order to determine the correct spec- 
tral dependences of the components of the microwave sky. 

Another possible explanation for the difficulty in fitting 
the data with our model is the presence of an unknown com- 
ponent. As already mentioned, several works have found an 
anomalous galactic emission at low frequency, which cor- 
relates with the thermal dust one. Although in principle a 
free-free component is fitting the data that we have used, its 
amplitude is higher than expected from estimations based on 
the Ha emission. The spectral index found for the free-free 
was also slightly higher than expected from theoretical pre- 
dictions. As a test, we have also run the code for a few mod- 
els with /3ff = —0.22, two of which gave slightly lower values 
of G than the chosen model. However, the synchrotron map 
was presenting negative features in the Galaxy, which is in- 
dicative of a degeneracy between the spectral parameters. 

To account for this anomalous component, we have 
tried to fit for a spinning dust component using the extra 
COBE-DIRBE channel in addition to the CMB, free-free, 
synchrotron and (thermal) dust emissions. The spectral pa- 
rameters have been fixed to the values found in §5.1. We have 
tried three different spinning dust models which correspond 
to the cold neutral medium, the warm neutral medium (de- 
scribed in Draine & Lazarian, 1998) and a combination of 
both. In the three considered cases, the G value was ap- 
preciably higher than in the optimal case and the CMB re- 
construction was lost. However, even if spinning dust were 
present in the data, it would have been very difficult to re- 
cover it. We are basically trying to reconstruct three com- 
ponents, CMB, free-free and spinning dust from just the 3 
COBE-DMR channels, which are very noisy, and therefore 
the data just do not have enough information. Moreover, a 
wider range of frequencies (between 5 and 60 GHz) would 
be necessary to distinguish between spinning dust and the 
other galactic components. 

Finally, we would like to point out that the combination 
of MEM and other reconstruction methods should also be 
investigated in the future. In particular, blind source sep- 
aration methods can infer the spectral dependence of the 
reconstructed components under certain assumptions. This 
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Figure 25. Top panel: CMB reconstructed map at 50 GHz 
smoothed with a 7 degree Gaussian beam for the case when 2 
dust components have been used. Middle panel: map resulting of 
subtracting the previous reconstruction from the one of Fig. 17. 
Bottom panel: total dust emission at 50 GHz obtained by adding 
the contribution from the two dust components smoothed with a 
7 degree Gaussian beam (non-linear scale). The CMB and dust 
reconstructions have been plotted in the same scale as the cor- 
responding reconstructions of Fig. 17. Units are (thermody- 
namical temperature). 

information could be used as input for the MEM algorithm, 
which would improve the results of each method alone. 
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APPENDIX A: 
DERIVATIVES 



FUNCTION AND 



In order to evaluate the x in real space we need to be able to 
predict the data at each pixel from the hidden variables h^m- 
The temperature field at a given pixel is normally written 
in terms of spherical harmonic coefficients. So, at a posi- 
tion i corresponding to angles {6, 4>) and at frequency v, the 
predicted noiseless data can be written as 



E 



E 



(Al) 



where aj"^ are the predicted harmonic coefficients at fre- 
quency v and the Ycm are the spherical harmonics: 

1/2 

P,„(cos0)e''"^. (A2) 



2£ + l (£-m)! 



47r {I + m) 



Taking into account that, for a real field, we have a(,-m = 
(— l)'"a|^, we can write dv as 

dl{6,(j>) = fcm W-jotm) cos(m0) 

1 = m=0 

-a (o^„) sin(m0) ] Ke.m (A3) 

where k.,n = 1 for rn = and km = 2 otherwise. A^m de- 
notes the factor that multiplies the exponential in the Yem 
expression. The are related to the hidden variables h^m 
by: 



(A4) 



where we have assumed that the beam is spherically sym- 
metric. B" is the beam in harmonic space at frequency i^, 
F is the conversion Uf x tic matrix, which encodes the fre- 
quency dependence of the components, is the Cholesky 
decomposition (an ric x ric matrix) and h^^ is the column 
vector containing the hidden modes for each of the mi- 
crowave components. The function is then written as 



(A5) 



2 [dg,,-<,,(h,^)]^ 

^ =2^2^ ^. ' 

U=l 1 = 1 ' 

where A^p is the number of pixels of each data map, d° ^ is 
the observed data at frequency f and pixel i and (T,/,i is the 
noise dispersion at frequency v and pixel i. 

If we denote and as the real and imaginary 
parts of the hidden harmonic coefficients for component c, 
we can write the first derivatives of as 



j9x! 
9rL 



"f d" -d" 

= -2^"^EE%^ 

,.=1 »=i 

d" -Jf 

= 2fc"^EE^^ 

u=i i=i 

X sm{m(t)i)KtmBi {F L)^c. 



(A6) 



The second derivatives are given by 



£l2 2 



2fe^ ^ ^ [cos{m(i>i)A.trnB'lf 



d2 2 
X 



1/=1 i=l 
x{FL)^,{FL),,,, 



^ ^ [sm(mcf)i)AimBef 



i/=l i=l ' 

x{FL),,{FL),, 



(A7) 



These second derivatives of x arc the main contribution to 
the second derivatives oi ip — x^ ^ <^S, since the term coming 
from the entropy is generally much smaller. They arc also 
the main contribution to the Hessian, since we are assum- 
ing that we can neglect the rest of the elements. Note that 
these derivatives are all positive and therefore (p will have 
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a well defined minimum in the space. This means that 
the Newton- Raphson (NR) method will always be driven in 
the right direction and should be able to find the desired 
global minimum. Note also that the second derivatives of 
depend only on the characteristics of the experiment, on the 
assumed spectral parameters and also in the initial power 
spectra but they are independent of the hem- Therefore we 
do not need to reevaluate these derivatives at each NR iter- 
ation. 
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